###########################################
# The Puzzle of Militia Containment (ISQ) #
# Figure Reproduction Code                #
# Brandon Bolte                           #
# 12/11/20                                #
###########################################

library(sandwich)
library(multiwayvcov)
library(lmtest)
library(MASS)
library(jtools)
library(ggplot2)
library(ggpubr)


setwd("C://Users/brand/Desktop/PGM_ISQ_Final_Files/Replication_Materials")
data <- read.csv("PGM_ISQmain.csv")

### Main models ###
logit11<-with(data=data, glm(cnew1plus~logtroop+rpr, family="binomial"))
vcov11<-cluster.vcov(logit11, data$war_id)
rob.fit11<-coeftest(logit11, vcov11, robust.se=T)
rob.fit11

logit12<-with(data=data, glm(cnew2plus~logtroop+rpr, family="binomial"))
vcov12<-cluster.vcov(logit12, data$war_id)
rob.fit12<-coeftest(logit12, vcov12, robust.se=T)
rob.fit12

logit21<-with(data=data, glm(cnew1plus~logtroop+rpr+shared_ethnic+criminals+peasants+semiofficial+jointmembership, family="binomial"))
vcov21<-cluster.vcov(logit21, data$war_id)
rob.fit21<-coeftest(logit21, vcov21, robust.se=T)
rob.fit21

logit22<-with(data=data, glm(cnew2plus~logtroop+rpr+shared_ethnic+criminals+peasants+semiofficial+jointmembership, family="binomial"))
vcov22<-cluster.vcov(logit22, data$war_id)
rob.fit22<-coeftest(logit22, vcov22, robust.se=T)
rob.fit22

logit31<-with(data=data, glm(cnew1plus~logtroop+rpr+shared_ethnic+criminals+peasants+semiofficial+jointmembership+contigwar+lngdpcapnew+separatist+polity+num_militias3, family="binomial"))
vcov31<-cluster.vcov(logit31, data$war_id)
rob.fit31<-coeftest(logit31, vcov31, robust.se=T)
rob.fit31

logit32<-with(data=data, glm(cnew2plus~logtroop+rpr+shared_ethnic+criminals+peasants+semiofficial+jointmembership+contigwar+lngdpcapnew+separatist+polity+num_militias3, family="binomial"))
vcov32<-cluster.vcov(logit32, data$war_id)
rob.fit32<-coeftest(logit32, vcov32, robust.se=T)
rob.fit32


#### FIGURES ####


### Main Manuscript ###

## Predicted Probability Plots ##

# Figure 1 - Logtroop #

troop1<-effect_plot(logit11, pred=logtroop,interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>1)",x.label="Weak Rebels (Troop Ratio)")+
  ylim(0,1)+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"))

troop2<-effect_plot(logit12, pred=logtroop, interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>2)",x.label="Weak Rebels (Troop Ratio)")+
  ylim(0,1)+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"))

ggarrange(troop1, troop2,
          labels=c("A", "B"),
          ncol=2, nrow=1)

ggsave("fig1.jpeg", 
       width=6.5, height=3, dpi=1200)



# Figure 2- RPR #

rpr1<-effect_plot(logit11, pred=rpr,interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>1)",x.label="Relative Political Reach")+
  ylim(0,1)+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"))

rpr2<-effect_plot(logit12, pred=rpr, interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>2)",x.label="Relative Political Reach")+
  ylim(0,1)+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12,face="bold"))


ggarrange(rpr1, rpr2,
          labels=c("A", "B"),
          ncol=2, nrow=1)

ggsave("fig2.jpeg", 
       width=6.5, height=3, dpi=1200)







### Appendix ###

## cox interaction plot (figure 1) ##

containment<-c(1,2,3)
rpreffects<-c(-0.679, -2.314, -3.949)
rprse<-c(0.833, 1.095, 2.019)
interaction1<-data.frame(cbind(containment, rpreffects, rprse))
interaction1$upper<-interaction1$rpreffects+1.65*interaction1$rprse
interaction1$lower<-interaction1$rpreffects-1.65*interaction1$rprse


ggplot(interaction1, aes(x=containment, y=rpreffects)) + geom_point(size=3) + 
  geom_errorbar(aes(ymin=lower, ymax=upper))+
  geom_hline(yintercept=0, linetype="dashed", color="darkgray", size=1)+
  labs(x="Containment", y="Effect of RPR (90% CIs)") +
  theme(text = element_text(size = 20))





## Predicted Probability Plots ##

#Model 3- LogTroop 
pdf("Model3_Rebels.pdf", width=10,height=6,paper='special')
effect_plot(logit21, pred=logtroop,interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>1)",x.label="Weak Rebels (Troop Ratio)")+
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()
#Model 3- RPR 
pdf("Model3_RPR.pdf", width=10,height=6,paper='special')

effect_plot(logit21, pred=rpr,interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>1)",x.label="Relative Political Reach")+
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()
#Model 4- LogTroop
pdf("Model4_Rebels.pdf", width=10,height=6,paper='special')

effect_plot(logit22, pred=logtroop, interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>2)",x.label="Weak Rebels (Troop Ratio)")+
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()
#Model 4- RPR
pdf("Model4_RPR.pdf", width=10,height=6,paper='special')

effect_plot(logit22, pred=rpr, interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>2)",x.label="Relative Political Reach")+
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()

#Model 5- LogTroop 
pdf("Model5_Rebels.pdf", width=10,height=6,paper='special')

effect_plot(logit31, pred=logtroop,interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>1)",x.label="Weak Rebels (Troop Ratio)")+
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()
#Model 5- RPR 
pdf("Model5_RPR.pdf", width=10,height=6,paper='special')

effect_plot(logit31, pred=rpr,interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>1)",x.label="Relative Political Reach")+
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()
#Model 6- LogTroop
pdf("Model6_Rebels.pdf", width=10,height=6,paper='special')

effect_plot(logit32, pred=logtroop, interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>2)",x.label="Weak Rebels (Troop Ratio)")+
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()
#Model 6- RPR
pdf("Model6_RPR.pdf", width=10,height=6,paper='special')

effect_plot(logit32, pred=rpr, interval=T, rug=T,rug.sides="b", robust=T, cluster="war_id",y.label="Pr(Containment>2)",x.label="Relative Political Reach")+
  
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=22,face="bold"))
dev.off()





